# ==============================================================================
# constants
# ==============================================================================

SPECLETTERS=['S','P','D','F','G','H','I','K','L','M','N','O','Q','R','T','U',\
             'V','W','X','Y','Z']

ATOMSYMBOLS=['A','H','He','Li','Be','B','C','N','O','F','Ne','Na','Mg','Al',\
             'Si','P','S','Cl','Ar','K','Ca','Sc','Ti','V','Cr','Mn','Fe','Co',\
             'Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y','Zr',\
             'Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I',\
             'Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy',\
             'Ho','Er','Tm','Yb','Lu','Hf','Ta','W','Re','Os','Ir','Pt','Au',\
             'Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa','U',\
             'Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No','Lr','Rf','Db',\
             'Sg','Bh','Hs','Mt','Ds','Rg','Uub','Uut','Uuq','Uup','Uuh','Uus',\
             'Uuo']

SUBSHELLS=['0s',\
           '1s',\
           '2s','2p',\
           '3s','3p','3d',\
           '4s','4p','4d','4f',\
           '5s','5p','5d','5f','5g',\
           '6s','6p','6d','6f','6g','6h',\
           '7s','7p','7d','7f','7g','7h','7i']

SSPIECES=[(0,0),\
          (1,0),\
          (2,0),(2,1),\
          (3,0),(3,1),(3,2),\
          (4,0),(4,1),(4,2),(4,3),\
          (5,0),(5,1),(5,2),(5,3),(5,4),\
          (6,0),(6,1),(6,2),(6,3),(6,4),(6,5),\
          (7,0),(7,1),(7,2),(7,3),(7,4),(7,5),(7,6)]

NSHELLS=['A','K','L','M','N','O','P','Q','R','S','T']

ROMAN=['0','I','II','III','IV','V','VI','VII','VIII','IX','X','XI','XII',\
       'XIII','XIV','XV','XVI','XVII','XVIII','XIX','XX','XXI','XXII',\
       'XXIII','XXIV','XXV','XXVI','XXVII','XXVIII','XXIX','XXX','XXXI',\
       'XXXII','XXXIII','XXXIV','XXXV','XXXVI','XXXVII','XXXVIII','XXXIX',\
       'XL']

EQUIV_LS_TERMS={}   # permitted LS terms for equivalent electrons; filled as 
                    # needed by rlf_equiv_terms()

RYDBERG=13.6056923
INVERSECM=109737.3
HC_RYDANG=911.26704
HC_EVNM=1239.8419
GF2A=0.66702

# ==============================================================================
# RAP exception type
# ==============================================================================

class ErrorRAP(Exception):
  def __init__(self,value):
    self.value=value
  def __str__(self):
    return repr(self.value)

# ==============================================================================
# RAP class helper functions
# ==============================================================================

def rlf_filetype(filename):
  """ 
  Check type of file: adf04, adf39l, TCCDW.DAT, olg, NIST, or xstar type 006.
  """
  if '.fits' in filename:
    rlf_have_pyfits(stop=True)
    try:
      a=pyfits.open(filename)
    except:
      pass
    else:
      origin=a[0].header['ORIGIN']
      dtype=a[1].header['HDUCLAS1']
      if origin == 'ATOMDB':
        if dtype == 'E_LEVEL':
          return 'atomdb_ev'
      a.close()

  fin=open(filename)
  lines=[]
  for line in fin:
    lines.append(line)
    if len(lines) == 10: break
  fin.close()

  if '=====' in lines[0]:     # must be TIPTOP
    if 'E(' in lines[1]: return 'tiptop'
    if 'WL(' in lines[1]: return 'tiptop_lines'

  for line in lines:
    if line[:3] == 'FAC':
      itype=int(lines[3].split()[-1])
      if itype == 1:
        return 'fac_lev'
      elif itype == 2:
        return 'fac_tr'
      else:
        break
    if line[:4] == 'SEQ=': return 'adf39l'
    tmp=line.split()
    if len(tmp) == 18 and tmp[-1] == '%' and int(tmp[0]) == 6: return 'x006'
    if len(tmp) in [14,15] and tmp[-1] == '%' and int(tmp[0]) == 50: 
      return 'x050'
    if len(tmp) == 12 and tmp[5] in SPECLETTERS: return 'chianti_elvlc'
    if len(tmp) == 10 and tmp[7] == '-': return 'chianti_wgfa'
    tmp=re.findall(r'(\d+)\s+(\d+)\s+(\d+)\s+(\d+\.)',line)
    if len(tmp) > 0:  return 'adf04'
    tmp=re.findall(r'(LS TERMS)',line)
    if len(tmp) > 0:  return 'tcc'
    if '**********' in line: return 'olg'
    if line.count('|') == 13: return 'nist_lines'
    if line.count('|') in [4,5]: return 'nist'
  raise ErrorRAP('cannot understand file type: '+filename)

# ------------------------------------------------------------------------------

def rlf_loadstates(filename,cup=None,**args):
  """ 
  Read file and create rap_states child object appropriate to detected or
  specified coupling scheme.  If coupling scheme is unspecified then try
  LSJ, LS, then CA.
  """
  if not os.path.exists(filename):
    msg='file, %s, not found'%filename
    raise ErrorRAP(msg)

  # coupling scheme is given...
  if cup == 'CA':
    try:
      out=rap_states_ca(filename,**args)
    except:
      msg='no CA state data found in %s'%filename
    else:
      return out
  if cup == 'LS':
    try:
      out=rap_states_ls(filename,**args)
    except:
      msg='no LS state data found in %s'%filename
    else:
      return out
  if cup == 'LSJ':
    try:
      out=rap_states_lsj(filename,**args)
    except:
      msg='no LSJ state data found in %s'%filename
    else:
      return out
  if cup is not None:
    msg='coupling scheme, %s, not understood'%cup
    raise ErrorRAP(msg)

  # coupling scheme is not given...
  try:
    out=rap_states_lsj(filename,**args)
  except ErrorRAP, e:
    msglsj=e.value
  else:
    return out
  try:
    out=rap_states_ls(filename,**args)
  except ErrorRAP, e:
    msgls=e.value
  else:
    return out
  try:
    out=rap_states_ca(filename,**args)
  except ErrorRAP, e:
    msgca=e.value
  else:
    return out
  msg=[]
  msg.append('LSJ: %s'%msglsj)
  msg.append('LS:  %s'%msgls)
  msg.append('CA:  %s'%msgca)
  msg.append('no state data found of any coupling scheme')
  msg='\n'.join(msg)
  raise ErrorRAP(msg)

# ==============================================================================
# functions acting on ion strings
# ==============================================================================

def rlf_iondecode(ion):
  """ 
  Return tuple, (z,ne), from ion string where z is the nuclear charge and ne is
  the number of electrons.
  """
  # look for sequence notation
  if '-like' in ion:
    try:
      tmp=ion.split('-like ')
      ne=ATOMSYMBOLS.index(tmp[0].strip().capitalize())
      z=ATOMSYMBOLS.index(tmp[1].strip().capitalize())
    except:
      msg="invalid ion string -- sequence notation"
      raise ErrorRAP(msg)
    return (z,ne)

  # look for astrophysical notation
  if ' ' in ion:
    tmp=ion.split()
    try:
      z=ATOMSYMBOLS.index(tmp[0].strip().capitalize())
      q=ROMAN.index(tmp[1].strip())-1
    except:
      msg="invalid ion string -- astrophysical notation"
      raise ErrorRAP(msg)
    else:
      return (z,z-q)
  # look for chemical notation
  tmp=re.findall(r'([A-Za-z][A-Za-z]?)\+?(\d+)\+?',ion)
  if not tmp:
    msg="invalid ion string -- chemical notation"
    raise ErrorRAP(msg)
  z=ATOMSYMBOLS.index(tmp[0][0].capitalize())
  q=int(tmp[0][1])
  return (z,z-q)

# ------------------------------------------------------------------------------

def rlf_ionencode(z,ne=None,form='chemical'):
  """ Take nuclear charge and electron count and return ion string. """
  if ne is None: (z,ne)=z
  q=z-ne
  x=ATOMSYMBOLS[z]
  if form[0].lower() == 'c':
    out=x+str(q)+'+'
  elif form[0].lower() == 'a':
    r=ROMAN[q+1]
    out=x+' '+r
  elif form[0].lower() == 's':
    out="%s-like %s"%(ATOMSYMBOLS[ne],ATOMSYMBOLS[z])
  else:
    msg="invalid output format"
    raise ErrorRAP(msg)
  return out

# ==============================================================================
# functions acting on sub-shells
# ==============================================================================

def rlf_subshell(i):
  """ Return sub-shell string of i'th sub-shell (in standard order). """
  if i < 1: 
    msg='sub-shell index out of range: i=%i'%i
    raise ErrorRAP(msg)
  if i <= len(SUBSHELLS)-1: return SUBSHELLS[i]
  n=int(round(sqrt(2.*float(i)-1.)))
  l=i-1-n*(n-1)/2
  if l > len(SPECLETTERS): 
    msg='angular momenta greater than '+str(len(SPECLETTERS)-1)+\
        ' not supported: i=%i => n,l=%i, %i'%(i,n,l)
    raise ErrorRAP(msg)
  return str(n)+SPECLETTERS[l].lower()

# ------------------------------------------------------------------------------

def rlf_ssindex(ss):
  """ 
  Return the standard order index of subshell, ss. 

  Supports continuum orbitals where the principal quantum number is replaced
  by 'k'.  The returned index is -(l+1) where l is the orbital angular 
  momentum.
  """
  if ss in SUBSHELLS: return SUBSHELLS.index(ss)
  if ss[0].upper() == 'K':                 # check for continuum orbital
    l=SPECLETTERS.index(ss[-1].upper())
    return -(l+1)
  tmp=re.findall(r'(\d+[A-Za-z])',ss)
  n=int(tmp[0][:-1])
  try:
    l=SPECLETTERS.index(tmp[0][-1].upper())
  except:
    msg='illegal symbol for L in spectroscopic notation: %s'%ss
    raise ErrorRAP(msg)
  if l > n-1:
    msg='illegal sub-shell: %s'%ss
    raise ErrorRAP(msg)
  i=n*(n-1)/2+l+1
  return i

# ------------------------------------------------------------------------------

def rlf_ssdecode(ssidx):
  """ Return tuple containing n,l for given sub-shell (string or index). """
  try:
    idx=int(ssidx)
  except:
    idx=rlf_ssindex(ssidx)
  try:
    out=SSPIECES[idx]
  except:
    n=int(round(sqrt(2.*float(idx)-1.)))
    l=idx-1-n*(n-1)/2
    out=(n,l)
  return out

# ------------------------------------------------------------------------------

def rlf_ssencode(n,l,index=False):
  """ Return sub-shell given n and l. """
  ss=str(n)+SPECLETTERS[l].lower()
  if index: ss=rlf_ssindex(ss)
  return ss

# ------------------------------------------------------------------------------

def rlf_ssoccupancy(ssidx):
  """ Return max occupancy of given sub-shell (string or index). """
  (n,l)=rlf_ssdecode(ssidx)
  return 4*l+2

# ------------------------------------------------------------------------------

def rlf_ss2jss(ssidx):
  """ Return list of allowed J subshells given a regular sub-shell. """
  (n,l)=rlf_ssdecode(ssidx)
  if l == 0: return [rlf_jssencode(n,l,+1)]
  return [rlf_jssencode(n,l,-1),rlf_jssencode(n,l,+1)]

# ------------------------------------------------------------------------------

def rlf_sssorter(x,y):
  """ 
  Use this function in the list sort() function in order to arrange items by
  their sub-shell order.
  """
  xx=rlf_ssindex(x)
  yy=rlf_ssindex(y)
  if xx < 0 and yy < 0:
    if xx < yy: return +1
    if xx > yy: return -1
    return 0
  if xx < 0: return +1
  if yy < 0: return -1
  if xx > yy: return +1
  if xx < yy: return -1
  return 0

# ------------------------------------------------------------------------------

def rlf_ss2eissner(idx):
  """ Return symbol for sub-shell in Eissner notation. """
  tidx=idx
  if not isinstance(tidx,int): tidx=rlf_ssindex(tidx)
  if idx <= 0: 
    msg='orbital not supported by Eissner notation; idx: %s'%idx
    raise ErrorRAP(msg)
  if tidx < 10:
    x=str(tidx)
  elif tidx < 36:
    x=chr(tidx-9-1+ord('A'))
  elif tidx < 62:
    x=chr(tidx-35-1+ord('a'))
  else:
    msg='orbital not supported by Eissner notation; idx: %s'%idx
    raise ErrorRAP(msg)
  return x

# ------------------------------------------------------------------------------

def rlf_eissner2ss(x,index=False):
  """ Convert sub-shell symbol from Eissner notation into standard form. """
  x=str(x)
  if len(x) > 1: 
    msg='Eissner orbital symbol too long: %s'%x
    raise ErrorRAP(msg)
  if x >= '1' and x <= '9':
    ss=int(x)
  elif x >= 'A' and x <= 'Z':
    ss=ord(x)-ord('A')+1+9
  elif x >= 'a' and x <= 'z':
    ss=ord(x)-ord('a')+1+35
  else:
    msg='invalid symbol for orbital in Eissner notation: %s'%x
    raise ErrorRAP(msg)
  if index: return ss
  return rlf_subshell(ss)

# ==============================================================================
# functions acting on JJ-coupled sub-shells
# ==============================================================================

def rlf_jsubshell(i):
  n=int(sqrt(float(i)-0.5))+1
  d=i-(n-1)**2
  l=d/2
  s=d%2
  out="%i%s"%(n,SPECLETTERS[l].lower())
  if s == 0:
    out+='-'
  else:
    out+='+'
  return out

# ------------------------------------------------------------------------------

def rlf_jssindex(ss):
  (n,l,s)=rlf_jssdecode(ss)
  return (n-1)**2+2*l+(s+1)/2

# ------------------------------------------------------------------------------

def rlf_jssdecode(iorss):
  if isinstance(iorss,int): iorss=rlf_jsubshell(iorss)
  try:
    n=int(iorss[:-2])
    l=SPECLETTERS.index(iorss[-2].upper())
    if '+' in iorss: 
      s=1
    elif '-' in iorss:
      s=-1
    else:
      raise
  except:
    msg='cannot understand J-subshell'
    raise ErrorRAP(msg)
  return (n,l,s)

# ------------------------------------------------------------------------------

def rlf_jssencode(n,l,s,index=False):
  if s not in [-1,+1]:
    msg='invalid value for sigma for J-subshell'
    raise ErrorRAP(msg)
  if l == 0 and s == -1:
    msg='sigma must be +1 when l=0'
    raise ErrorRAP(msg)
  idx=(n-1)**2+2*l+(s+1)/2
  if index: return idx
  return rlf_jsubshell(idx)

# ------------------------------------------------------------------------------

def rlf_jssoccupancy(iorss):
  (n,l,s)=rlf_jssdecode(iorss)
  return 2*l+s+1

# ------------------------------------------------------------------------------

def rlf_jsssorter(x,y):
  xx=rlf_jssindex(x)
  yy=rlf_jssindex(y)
  return cmp(xx,yy)

# ==============================================================================
# functions acting term strings
# ==============================================================================

def rlf_termsorter(x,y):
  """
  Use this function in the list sort() in order to arrange entries according 
  to term order, e.g. tlist.sort(rlf_termsorter).
  """
  return cmp(rlf_termdecode(x),rlf_termdecode(y))

# ------------------------------------------------------------------------------

def rlf_termdecode(term):
  """ Convert term string to tuple of 2S+1 and L. """
  try:
    g=int(term[:-1])
    l=SPECLETTERS.index(term[-1])
  except:
    msg='bad term string: %s'%term
    raise ErrorRAP(msg)
  return (g,l)

# ------------------------------------------------------------------------------

def rlf_termencode(g,l):
  """ Return term string given 2S+1 and L. """
  if l < 0 or l > len(SPECLETTERS)-1:
    msg='angular momentum out of supported range; L=%s'%l
    raise ErrorRAP(msg)
  return str(g)+SPECLETTERS[l]

# ------------------------------------------------------------------------------

def rlf_termweight(term):
  """ Return term weight: (2S+1)(2L+1). """
  (g,l)=rlf_termdecode(term)
  return g*(2*l+1)

# ------------------------------------------------------------------------------

def rlf_term_2jvals(term):
  """ Return list of allowed 2J values for given term. """
  (g,l)=rlf_termdecode(term)
  j2min=abs(g-(2*l+1))
  j2max=g+(2*l+1)-2
  return range(j2min,j2max+1,2)

# ------------------------------------------------------------------------------

def _rlf_calc_equiv_lsterms(lorb,ne):
  """ 
  Calculate dictionary of LS terms of orbital with angular momentum, lorb, 
  and having ne electrons.

  The returned dictionary has keys composed of tuples having the form 
  (2S+1,L,p), where 2S+1,L is the usual term identifiers and p is the parity.
  The value of the dictionary with a given key is the number of terms in that
  list with that key.
  """
  def to_binary(x,npad=0):
    """ Convert integer into binary representation via list """
    xx=int(x)
    out=[]
    while xx > 0:
      (xx,y)=divmod(xx,2)
      out.append(y)
    if npad > len(out):
      for i in range(len(out),npad): out.append(0)
    out.reverse()
    return out
  
  def array_mult(x,y):
    """ Multiply elements of two arrays together to form new array """
    if len(x) != len(y): raise ErrorRAP('arrays must be of equal size')
    out=[]
    for i in range(len(x)):
      out.append(x[i]*y[i])
    return out
  
  def dot_prod(x,y):
    """ Take dot-product of two arrays of equal size """
    vec=array_mult(x,y)
    return sum(vec)

  # short-cut for single electron (hole)
  if ne == 1 or ne == 4*lorb+1:
    termlist={}
    key=(2,lorb,(1-(-1)**lorb)/2)
    termlist[key]=1
    return termlist

  # define orbital
  ll=lorb
  nelec=ne
  parity=(1-(-1)**(ll*nelec))/2
  
  # form arrays defining angular momenta and spins values
  lval=[]
  for i in range(2*ll+1):
    lval.append(-ll+i)
  lval*=2
  sval=[-1]*(2*ll+1)+[+1]*(2*ll+1)
  
  # calculate all possible combinations of L and S
  nbox=4*ll+2
  nmax=2**nbox-1
  stot=[]
  ltot=[]
  for n in range(1,nmax+1):
    arr=to_binary(n,npad=nbox)
    if sum(arr) == nelec:
      stot.append(dot_prod(sval,arr))
      ltot.append(dot_prod(lval,arr))

  # form matrix with frequency distribution of L and S
  lmax=max(ltot)
  smax=max(stot)
  col=[]
  l=-lmax
  while l <= lmax:
    col.append(0)
    l+=1
  table=[]
  s=-smax
  while s <= smax:
    table.append(col[:])
    s+=2
  for i in range(len(ltot)):
    jl=ltot[i]+lmax
    js=(stot[i]+smax)/2
    table[js][jl]+=1

  # get list of terms
  termlist={}
  jlmax=(len(col)+1)/2-1
  jsmax=(len(table)+1)/2-1
  js=0
  while js <= jsmax:
    jl=0
    s=2*js-smax
    while jl <= jlmax:
      l=jl-lmax
      n=table[js][jl]
      if n > 0:
        key=(abs(s)+1,abs(l),parity)
        for kl in range(jl,len(col)-jl):
          for ks in range(js,len(table)-js):
            table[ks][kl]=table[ks][kl]-n
        termlist[key]=n
      jl=jl+1
    js=js+1
  return termlist

# ------------------------------------------------------------------------------

def _rlf_lsterm_couple(t1,t2,numterms=1):
  """ 
  Internal function for coupling two terms together; each term is a tuple of
  the form, (2S+1,L,p).  Returns a LS term list with the number of terms set
  to numterms.
  """
  (s1,l1,p1)=t1
  (s2,l2,p2)=t2
  n=int(numterms)
  smin=abs(s1-s2)
  smax=abs(s1+s2-2)
  lmin=abs(l1-l2)
  lmax=abs(l1+l2)
  p=p1+p2
  if p == 2: p=0
  termlist={}
  for s in range(smin,smax+1,2):
    for l in range(lmin,lmax+1):
      nt=(s+1,l,p)
      termlist[nt]=n
  return termlist

# ------------------------------------------------------------------------------

def rlf_equiv_lsterms(lorb,ne):
  """ 
  Return dictionary of LS terms from an orbital with angular momentum, lorb, 
  having ne equivalent electrons.

  The dictionary of terms is only calculated once and then stored in the 
  EQUIV_LS_TERMS global variable (reset after session ends).
  """
  try:
    ll=int(lorb)
    nn=int(ne)
  except:
    return False
  if nn > 4*ll+2: 
    msg='too many electrons for shell: lorb,ne: %s, %s'%(lorb,ne)
    raise ErrorRAP(msg)
  ss=(ll,nn)
  if ss in EQUIV_LS_TERMS.keys(): return EQUIV_LS_TERMS[ss]
  termlist=_rlf_calc_equiv_lsterms(ll,nn)
  EQUIV_LS_TERMS[ss]=termlist
  return termlist

# ------------------------------------------------------------------------------

def rlf_lsterm_merge(termlist1,termlist2):
  """ Merge termlists together and return result. """
  if not isinstance(termlist1,dict):
    msg='input term lists must be dictionaries, received: %s'%type(termlist1)
    raise ErrorRAP(msg)
  if not isinstance(termlist2,dict):
    msg='input term lists must be dictionaries, received: %s'%type(termlist2)
    raise ErrorRAP(msg)
  termlist=termlist1.copy()
  for t,n in termlist2.iteritems():
    if t not in termlist:
      termlist[t]=n
    else:
      termlist[t]+=n
  return termlist

# ------------------------------------------------------------------------------

def rlf_lsterm_couple(termlist1,termlist2):
  """ 
  Couple all terms from termlist1 with terms from termlist2 and return the
  resulting dictionary.
  """
  if not isinstance(termlist1,dict):
    msg='input term lists must be dictionaries, received: %s'%type(termlist1)
    raise ErrorRAP(msg)
  if not isinstance(termlist2,dict):
    msg='input term lists must be dictionaries, received: %s'%type(termlist2)
    raise ErrorRAP(msg)
  termlist={}
  for t1,n1 in termlist1.iteritems():
    for t2,n2 in termlist2.iteritems():
      n=n1*n2
      termlist=rlf_lsterm_merge(termlist,_rlf_lsterm_couple(t1,t2,n))
  return termlist

# ------------------------------------------------------------------------------

def rlf_num_lsterms(termlist):
  """ Return number of terms in list """
  if not isinstance(termlist,dict):
    msg='input term list must be a dictionary, received: %s'%type(termlist)
    raise ErrorRAP(msg)
  total=0
  for n in termlist.itervalues(): 
    total+=n
  return total

# ==============================================================================
# general configuration functions
# ==============================================================================

def rlf_confsorter(x,y):
  """
  Use this function in the list sort() in order to arrange entries according 
  to configuration order, e.g. alist.sort(rlf_confsorter).
  """
  if not isinstance(x,rap_config):
    msg='only rap_config types can be sorted; received: %s'%type(x)
    raise ErrorRAP(msg)
  if not isinstance(y,rap_config):
    msg='only rap_config types can be sorted; received: %s'%type(y)
    raise ErrorRAP(msg)
  ssx=x.subshells(indices=True)
  ssy=y.subshells(indices=True)
  if x.num_electrons() < y.num_electrons(): return -1
  if x.num_electrons() > y.num_electrons(): return +1
  ssmax=max(max(ssx),max(ssy))
  for i in range(1,ssmax+1):
    if i not in ssx and i not in ssy: continue
    if i not in ssx: return +1
    if i not in ssy: return -1
    wx=x.occupancy(i)
    wy=y.occupancy(i)
    if wx < wy: return +1
    if wx > wy: return -1
    if i == ssmax: return 0
  raise ErrorRAP('sanity check - should have returned a value')

# ------------------------------------------------------------------------------

def rlf_jconfsorter(x,y):
  """
  Use this function in the list sort() in order to arrange entries according 
  to configuration order, e.g. alist.sort(rlf_confsorter).
  """
  if not isinstance(x,rap_jconfig):
    msg='only rap_jconfig types can be sorted; received: %s'%type(x)
    raise ErrorRAP(msg)
  if not isinstance(y,rap_jconfig):
    msg='only rap_jconfig types can be sorted; received: %s'%type(y)
    raise ErrorRAP(msg)
  ssx=x.subshells(indices=True)
  ssy=y.subshells(indices=True)
  if x.num_electrons() < y.num_electrons(): return -1
  if x.num_electrons() > y.num_electrons(): return +1
  ssmax=max(max(ssx),max(ssy))
  for i in range(1,ssmax+1):
    if i not in ssx and i not in ssy: continue
    if i not in ssx: return +1
    if i not in ssy: return -1
    wx=x.occupancy(i)
    wy=y.occupancy(i)
    if wx < wy: return +1
    if wx > wy: return -1
    if i == ssmax: return 0
  raise ErrorRAP('sanity check - should have returned a value')

# ------------------------------------------------------------------------------

def rlf_genconf(aconf,force=False):
  """
  Generate list of configurations from a configuration string which includes
  n- and l- variables.  

  Example: '2l5 nl n=3,4 l=0-2' will expand into a 
  list of 18 configurations.  Multiple variables can be used as long as the
  n-variables begin with 'n' and the l-variables begin with 'l'.  Also, the
  variables cannot contain any numbers.  Finally, to avoid confusion, the 
  n-variable cannot contain any l's.  For each n-variable, a range must be
  specified (as in the above example).  The l-variables will take all values
  by default unless a range is given.  For cases where an l-electron is
  permissible (L=8 electrons), the routine will assume you want the single
  electron unless the force variable is True or 'l' is used as a variable 
  somewhere else in the string.  Note that no spaces are allowed in the 
  range definitions (c.f. n = 3,4 is legal; n=3, 4 is illegal).  

  Examples:  
    2l5                          ( 3 configurations)
    3s 3l                        ( 3 configurations)
    nl npl n=3,4 np=3-5 l=0-2    (54 configurations)
  """
  def fill_conf(n,lvals,ne,conf0=''):
  # A recursive function to generate a list of configurations with ne 
  # electrons for a given n-value for values of l in lvals.  The variable
  # conf0 is used in the recursive calling process.
    ll=lvals[0]
    if len(lvals) == 1:
      if ne > 4*ll+2: return None
      if SPECLETTERS[ll].lower() == 'l':
        ss=str(n)+'LTEMP'+str(ne)
      else:
        ss=str(n)+SPECLETTERS[ll].lower()+str(ne)
      return [(conf0+' '+ss).strip()]
    else:
      out=[]
      for k in range(ne+1):
        tne=ne-k
        if tne > 4*ll+2: continue
        if SPECLETTERS[ll].lower() == 'l':
          ss=str(n)+'LTEMP'+str(tne)
        else:
          ss=str(n)+SPECLETTERS[ll].lower()+str(tne)
        tconf=(conf0+' '+ss).strip()
        tmp=fill_conf(n,lvals[1:],k,tconf)
        if isinstance(tmp,list):
          out.extend(tmp)
        elif tmp:
          out.append(tmp.strip())
      return out

  def parse_range(x):
  # Return a list of values from the user input range of n or l.
  # (Example: 3,5-7 -> [3,5,6,7] )
    tmp1=re.findall(r'^\d+',x)
    tmp2=re.findall(r'([,-])(\d+)',x)
    if len(tmp1) == 0: return False
    s=[tmp1[0]]
    gotdash=False
    if len(tmp2) != 0:
      for t in tmp2:
        if t[0] == ',':
          s.append(t[1])
          gotdash=False
        else:
          if gotdash == True: return False
          s[-1]=str(s[-1])+t[0]+t[1]
          gotdash=True
    out=[]
    for it in s:
      if '-' not in it:
        out.append(int(it))
      else:
        tmp=re.findall(r'(\d+)',it)
        out.extend(range(int(tmp[0]),int(tmp[1])+1))
    return out
  
  # count number of subshells
  tmp=re.findall(r'([^\s]+)',str(aconf))
  sscount=0
  for i in range(len(tmp)):
    t=tmp[i]
    if '=' not in t: 
      sscount=sscount+1
      try:
        ne=int(t[-1])
      except:
        tmp[i]=t+'1'
  conf=''
  for t in tmp: conf=conf+' '+t
  conf=conf.strip()

  # identify metavariables
  tnmeta=[]
  tlmeta=[]
  ce=re.findall(r'(n[A-Za-km-z]*|\d+)(l[A-Za-z]*|[a-z])(\d+)',conf)
  if len(ce) == 0: 
    msg='invalid configuration: %s'%conf
    raise ErrorRAP(msg)
  for ss in ce:
    # check if 'n' is an integer
    try:
      n=int(ss[0])
    except:
      if ss[0] not in tnmeta: tnmeta.append(ss[0])
      n=0
    # check if 'l' is spectroscopic letter
    if ss[1][0] != 'l':
      if ss[1].upper() not in SPECLETTERS: 
        msg="l metavariable must start with 'l': %s"%ss
        raise ErrorRAP(msg)
      l=SPECLETTERS.index(ss[1].upper())
    else:
      l=-1
      if len(ss[1]) > 1:
        if ss[1] not in tlmeta: tlmeta.append(ss[1])
      # if l='l', check if it is legal for given n
      # always metavar if force is True
      else:
        if n < 9 or force == True: 
          if ss[1] not in tlmeta: tlmeta.append(ss[1])
    # check if legal number of electrons
    ne=int(ss[2])
    if l >= 0:
      if ne > 4*l+2: 
        msg='too many electrons for l=%i'%l
        raise ErrorRAP(msg)
    elif n > 0:
      if ne > 2*n*n: 
        msg='too many electrons for n=%i'%n
        raise ErrorRAP(msg)
  # check that all subshells were legal
  if len(ce) != sscount:
    msg='invalid configuration string: %s'%conf
    raise ErrorRAP(msg)

  # parse range information
  nmeta={}
  lmeta={}
  tmp=re.findall(r'([^\s]+)=([^\s]+)',conf)
  for t in tmp:
    if t[0] in tnmeta:
      r=parse_range(t[1])
      if r != False: nmeta[t[0]]=r
    elif t[0] in tlmeta:
      r=parse_range(t[1])
      if r != False: lmeta[t[0]]=r
    else:
      continue
  for n in tnmeta:
    if n not in nmeta:
      msg='must specify range for all n metavariables'
      raise ErrorRAP(msg)
  for l in tlmeta:
    if l not in lmeta:
      lmeta[l]=[]

  # get list of possibilities for each subshell
  nmin=None
  lmin=None
  poss=[]
  null=[]
  for ss in ce:
    if ss[0] in nmeta:
      nval=nmeta[ss[0]]
    else:
      nval=[int(ss[0])]
    if nmin is None: nmin=min(nval)
    if ss[1] in lmeta:
      lval=lmeta[ss[1]]
    else:
      lval=[SPECLETTERS.index(ss[1].upper())]
    tlmin=min(lval,0)
    if lmin is None: lmin=tlmin
    if min(nval) == nmin:
      if tlmin < lmin: lmin=tlmin
    elif min(nval) < nmin:
      lmin=tlmin
      nmin=min(nval)
    ne=ss[2]
    p=[]
    for n in nval:
      if len(lval) == 0:
        lrange=range(int(n))
      else:
        lrange=[]
        for la in lval:
          if la < int(n): lrange.append(la)
      if(int(ne)) == 0:
        for l in lrange: null.append(str(n)+SPECLETTERS[l].lower())
      p.extend(fill_conf(int(n),lrange,int(ne)))
    poss.append(p)

  # remove temporary variable LTEMP
  for i in range(len(poss)):
    for j in range(len(poss[i])):
      try:
        poss[i][j]=poss[i][j].replace('LTEMP','l')
      except:
        pass

  # construct configurations
  out=poss.pop(0)
  for p in poss:
    tmp=[]
    for ap in p:
      for ss in out:
        tmp.append(ss+' '+ap)
    out=tmp

  # form and sort configuration
  nout=[]
  ne=None
  for i in range(len(out)): 
    try:
      c=rap_config(out[i],ne=ne)
    except ErrorRAP, e:
      if 'maximum occupancy exceeded' in e.value: continue
    nout.append(c)
    if ne is None: ne=nout[-1].num_electrons()
  nout.sort(rlf_confsorter)
  return nout

# ==============================================================================
# angular algebra
# ==============================================================================

def pyang_triangle(j1,j2,j3):
  """ Return True or False if triangle relations are met. """
  if j1+j2 < j3: return False
  if j2+j3 < j1: return False
  if j3+j1 < j2: return False
  return True

# ------------------------------------------------------------------------------

def pyang_wigner3j(j1,j2,j3,m1,m2,m3):
  """ Return Wigner 3-j symbol. """
  if m1+m2+m3 != 0.0: return 0.0
  if int(2*(j1+j2+j3)) % 2 != 0: return 0.0
  if not pyang_triangle(j1,j2,j3): return 0.0

  out=factorial(j1+j2-j3)*factorial(j1-j2+j3)*factorial(-j1+j2+j3)
  out*=factorial(j1-m1)*factorial(j1+m1)
  out*=factorial(j2-m2)*factorial(j2+m2)
  out*=factorial(j3-m3)*factorial(j3+m3)
  out=float(out)/float(factorial(j1+j2+j3+1))
  out=sqrt(out)

  kmin=int(max(0,j2-j3-m1,j1-j3+m2))
  kmax=int(min(j1+j2-j3,j1-m1,j2+m2))
  ksum=0.0
  for k in range(kmin,kmax+1):
    t=factorial(k)*factorial(j1+j2-j3-k)*factorial(j1-m1-k)
    t*=factorial(j2+m2-k)*factorial(j3-j2+m1+k)*factorial(j3-j1-m2+k)
    ksum+=(-1.)**k/float(t)

  out*=(-1.)**int(j1-j2-m3)*ksum
  return out

# ------------------------------------------------------------------------------

def pyang_wigner6j(j1,j2,j3,l1,l2,l3):
  """ Return Wigner 6-j symbol. """

  def tri(a,b,c):
    if int(2*(a+b+c)) % 2 != 0: return 0.0
    if not pyang_triangle(a,b,c): return 0.0
    out=factorial(a+b-c)*factorial(a-b+c)*factorial(-a+b+c)
    out=float(out)/float(factorial(a+b+c+1.))
    out=sqrt(out)
    return out

  out=tri(j1,j2,j3)*tri(j1,l2,l3)*tri(l1,j2,l3)*tri(l1,l2,j3)
  if out == 0.0: return out

  kmin=int(max(j1+j2+j3,j1+l2+l3,l1+j2+l3,l1+l2+j3))
  kmax=int(min(j1+j2+l1+l2,j2+j3+l2+l3,j3+j1+l3+l1))
  ksum=0.0
  for k in range(kmin,kmax+1):
    t=factorial(k-j1-j2-j3)*factorial(k-j1-l2-l3)*factorial(k-l1-j2-l3)
    t*=factorial(k-l1-l2-j3)*factorial(j1+j2+l1+l2-k)
    t*=factorial(j2+j3+l2+l3-k)*factorial(j3+j1+l3+l1-k)
    ksum+=(-1.)**k*float(factorial(k+1))/float(t)

  out*=ksum
  return out

# ------------------------------------------------------------------------------

def pyang_wigner9j(j1,j2,j3,j4,j5,j6,j7,j8,j9):
  """ Return Wigner 9-j symbol. """
  j2min=int(2.*max(abs(j1-j9),abs(j4-j8),abs(j2-j6)))
  j2max=int(2.*max(j1+j9,j4+j8,j2+j6))
  jsum=0.0
  for jj2 in range(j2min,j2max+1,2):
    j=0.5*float(jj2)
    t=(-1.)**(jj2)*float(jj2+1)*pyang_wigner6j(j1,j4,j7,j8,j9,j)
    t*=pyang_wigner6j(j2,j5,j8,j4,j,j6)*pyang_wigner6j(j3,j6,j9,j,j1,j2)
    jsum+=t
  return jsum

# ==============================================================================
# numerical functions
# ==============================================================================

def rlf_aitken(x,y,x0,logx=False,logy=False):
    """
    Use aitken method to interpolate data points, (x_i,y_i) in input arrays x 
    and y at point x0.
    """
    if not isinstance(x,list) and not isinstance(x,tuple):
      msg='input arrays must be lists; received: %s'%type(x)
      raise ErrorRAP(msg)
    if not isinstance(y,list) and not isinstance(y,tuple):
      msg='input arrays must be lists; received: %s'%type(y)
      raise ErrorRAP(msg)
    if len(x) != len(y):
      msg='input arrays have different lengths: len(x),len(y)=%s,%s'% \
          (len(x),len(y))
      raise ErrorRAP(msg)
    if logx:
        tx=[]
        for xx in x: tx.append(log10(xx))
        tx0=log10(x0)
    else:
        tx=x[:]
        tx0=x0
    if logy:
        ty=[]
        for yy in y: ty.append(log10(yy))
    else:
        ty=y[:]
    npts=len(tx)
    h=[]
    yt=[]
    for k in range(npts):
        h.append(tx0-tx[k])
        yt.append(ty[k])
    for k in range(npts-1):
        for l in range(k+1,npts):
            yt[l]=(h[k]*yt[l]-h[l]*yt[k])/(tx[l]-tx[k])
    out=yt[npts-1]
    if logy: out=10.0**out
    return out

# ==============================================================================
# xstar functions
# ==============================================================================

def xstar_ionencode(zeff,ne):
  """ Return xstar ion number. """
  if zeff < 1:
    msg='nuclear charge not physical: Z=%s'%zeff
    raise ErrorRAP(msg)
  if ne > zeff:
    msg='negative ions not supported: Z,ne:%s,%s'%(zeff,ne)
    raise ErrorRAP(msg)
  return zeff*(zeff+1)/2+1-ne

# ------------------------------------------------------------------------------

def xstar_iondecode(ionno):
  """ Return tuple, (Z,ne), for the given xstar ion number. """
  if ionno < 1:
    msg='xstar ion number must be positive: ionno=%s'%ionno
    raise ErrorRAP(msg)
  zeff=int(ceil(-0.5+0.5*sqrt(8.*float(ionno)-3.)))
  ne=zeff*(zeff+1)/2+1-ionno
  return (zeff,ne)

# ==============================================================================
# ADAS functions
# ==============================================================================

def rlf_adas2float(x):
  """ 
  Takes a number in an ADAS file (no 'E' given) and gives a floating point 
  number.
  """
  return float(x[:-3]+'E'+x[-3:])

# ==============================================================================
# functions acting AUTOSTRUCTURE olg files
# ==============================================================================

def rlf_olg_params(filename,indices=False,terse=False):
  """
  Return a dictionary of Thomas-Fermi scaling parameters in the given 
  AUTOSTRUCTURE output file where the key is the subshell indentifier, eg. 2p.
  If terse then return a list of scaling parameters with no indentifier.
  """
  if rlf_filetype(filename) != 'olg': 
    msg="invalid file type for %s"%filename
    raise ErrorRAP(msg)
  fin=open(filename,'r')
  while 1:
    line=fin.readline()
    if line == '': 
      msg="end-of-file reached before orbital specification"
      raise ErrorRAP(msg)
    if 'INTERNAL ORBITALS' in line: break
  torbs=[]
  while 1:
    line=fin.readline()
    if 'CONFIGURATION' in line: break
    if 'N L' in line: torbs.extend(line.split()[2:])
  norbs=len(torbs)/2
  orbs=[]
  for i in range(norbs):
    ss=str(torbs[2*i])+SPECLETTERS[int(torbs[2*i+1])].lower()
    if indices: ss=rlf_ssindex(ss)
    orbs.append(ss)
  while 1:
    line=fin.readline()
    if line == '':
      msg="end-of-file reached before scaling parameters"
      raise ErrorRAP(msg)
    if 'SCALING' in line: break
  bline=line[:-1]
  while 1:
    line=fin.readline()
    if line == '' or 'I' in line: break
    bline+=line[:-1]
  tmp=re.findall(r'(\d+\.\d+)',bline)
  params=[]
  for t in tmp: params.append(float(t))
  fin.close()
  if terse: return params
  if len(orbs) != len(params):
    msg="number of scaling parameters does not match number of orbitals"
    raise ErrorRAP(msg)
  return dict(zip(orbs,params))

# ==============================================================================
# compatibility functions
# ==============================================================================

def rlf_have_pyfits(stop=False):
  """ Return True if the pyfits module is installed. """
  if 'pyfits' in sys.modules.keys(): return True
  if stop:
    msg='module, pyfits, is required'
    raise ErrorRAP(msg)
  return False

# ------------------------------------------------------------------------------

# implement sorted() function for Python v2.3
try:
  sorted
except NameError: # doesn't have "sorted"
  def sorted(arr):
    out=arr[:]
    out.sort()
    return out

# ------------------------------------------------------------------------------

# implement reversed() function for Python v2.3
try:
  reversed
except NameError: # doesn't have "reversed"
  def reversed(arr):
    out=arr[:]
    out.reverse()
    return out

# ------------------------------------------------------------------------------

# implement itertools.product
try:
  iter_product=itertools.product
except:
  def iter_product(*args, **kwds):
    # product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
    # product(range(2), repeat=3) --> 000 001 010 011 100 101 110 111
    pools = map(tuple, args) * kwds.get('repeat', 1)
    result = [[]]
    for pool in pools:
      result = [x+[y] for x in result for y in pool]
    for prod in result:
      yield tuple(prod)

# ------------------------------------------------------------------------------

# implement itertools.combinations
try:
  iter_combinations=itertools.combinations
except:
  def iter_combinations(iterable, r):
    # combinations('ABCD', 2) --> AB AC AD BC BD CD
    # combinations(range(4), 3) --> 012 013 023 123
    pool = tuple(iterable)
    n = len(pool)
    if r > n:
      return
    indices = range(r)
    yield tuple([pool[i] for i in indices])
    while True:
      for i in reversed(range(r)):
        if indices[i] != i + n - r:
          break
      else:
        return
      indices[i] += 1
      for j in range(i+1, r):
        indices[j] = indices[j-1] + 1
      yield tuple([pool[i] for i in indices])

# ------------------------------------------------------------------------------

# implement itertools.combinations_with_replacement
try:
  iter_combinations_with_replacement=itertools.combinations_with_replacement
except:
  def iter_combinations_with_replacement(iterable, r):
    """ Python 2.7 function in itertools. """
    pool = tuple(iterable)
    n = len(pool)
    for indices in iter_product(range(n), repeat=r):
      if sorted(indices) == list(indices):
        yield tuple([pool[i] for i in indices])



# ==============================================================================
# miscellaneous functions
# ==============================================================================

def rlf_termsize():
  """
  Return tuple containing the height and width of the terminal window.
  """
  return struct.unpack('hh',fcntl.ioctl(0, termios.TIOCGWINSZ, '1234'))

# ------------------------------------------------------------------------------

def rlf_equiv(alist):
  """ Return True if all items in list, alist, are identical. """
  if alist.count(alist[0]) == len(alist): return True
  return False

# ------------------------------------------------------------------------------

def rlf_binom(a,b):
  """ Return binomial of (a b). """
  if b > a: raise ErrorRAP('invalid binomial coefficients: b > a')
  tb=max(b,a-b)
  top=1.0
  for i in range(tb+1,a+1): top*=i
  bot=1.0
  for i in range(1,a-tb+1): bot*=i
  return top/bot

# ------------------------------------------------------------------------------

def rlf_ordinal(i):
  """ Return ordinal of given integer, e.g. 1 yields 1st. """
  if i%10 == 1 and i%100 != 11:
    return str(i)+'st'
  elif i%10 == 2 and i%100 != 12:
    return str(i)+'nd'
  elif i%10 == 3 and i%100 != 13:
    return str(i)+'rd'
  else:
    return str(i)+'th'

# ------------------------------------------------------------------------------

def rlf_wraptext(text,length,indent=0):
  """
  Take string and return a list of strings where the length of each is at most
  the given length.  Separation of the initial string will only occur at 
  spaces.  An optional indentation can be given for all lines following the
  first.
  """
  if len(text) < length: return text
  out=[]
  pieces=text.split()
  line=pieces.pop(0)
  while len(pieces) > 0:
    if len(line)+len(pieces[0]) > length:
      out.append(line)
      line=" "*indent+pieces.pop(0)
    else:
      line+=" "+pieces.pop(0)
  out.append(line)
  return out

